LinkedIn GitHub

  تمرین  17  شبیه  سازی استاد جباری   نیمسال جاری

این کد رو درست کن

بخش اول: روش تقریبی (Approximate Method)

مفهوم کلی: قضیه حد مرکزی

کل ایده این بخش این است که گاهی تولید دقیق اعداد تصادفی (مثلاً برای توزیع نرمال یا پواسون) سخت یا زمان‌بر است. در عوض، ما از قوانین ریاضی مانند قضیه حد مرکزی (CLT) استفاده می‌کنیم تا با روشی ساده‌تر و سریع‌تر، اعدادی تولید کنیم که تقریباً همان توزیع را دارند.

اگر تعداد زیادی نمونه تصادفی (\(n\) بزرگ) از هر توزیعی داشته باشیم، میانگین (یا مجموع) آن‌ها شبیه به توزیع نرمال رفتار می‌کند.

فرمول پایه و شبیه‌سازی نرمال استاندارد

اگر \(X_1, ..., X_n\) متغیرهای تصادفی با میانگین \(\mu\) و واریانس \(\sigma^2\) باشند، برای \(n\) بزرگ داریم:

\[Z = \frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}} \approx N(0,1)\]

تکنیک طلایی (\(n=12\)): برای تولید توزیع نرمال استاندارد (\(Z\))، ساده‌ترین راه استفاده از ۱۲ عدد تصادفی یکنواخت \(U(0,1)\) است. چرا ۱۲؟ چون واریانس توزیع یکنواخت برابر \(\frac{1}{12}\) است و محاسبات بسیار ساده می‌شود:

\[Z = \sum_{i=1}^{12} U_i - 6 \approx N(0,1)\]

شبیه‌سازی توزیع پواسون (\(P(\lambda)\))

برای تولید عدد تصادفی پواسون، از تقریب دوجمله‌ای به نرمال یا پواسون استفاده می‌کنیم. اگر در توزیع دوجمله‌ای \(B(n,p)\)، مقدار \(n\) بزرگ و \(p\) کوچک باشد، توزیع شبیه پواسون با \(\lambda = np\) می‌شود.

الگوریتم: ۱. محاسبه مقدار \(n\): فرمول \(n = \max(30, \lceil 20\lambda \rceil)\). ۲. تولید \(n\) عدد تصادفی یکنواخت. ۳. شمارش اعدادی که کمتر یا مساوی \(p=0.05\) هستند.

فرمول ریاضی شمارش: \[x = \sum_{i=1}^{n} I(u_i \le 0.05)\]

پیاده‌سازی در R

در کد زیر، تابع sim_poisson دقیقاً بر اساس الگوریتم تقریبی فوق پیاده‌سازی شده است:

sim_poisson <- function(L,p) {
  
  
  p <- p
  
  n <- max(30, (20*L))
  
  
  X <- rbinom(n = 1, size = n, prob = p)
  
  return(X)
}

replicate(10, sim_poisson(L = 3,0.05))
##  [1] 4 3 4 3 5 3 1 5 1 3

بخش دوم: روش رد و پذیرش (Acceptance-Rejection)

معرفی روش

این روش یکی از تکنیک‌های اصلی در شبیه‌سازی آماری برای تولید اعداد تصادفی از توزیع‌های پیچیده است. زمانی که نمی‌توانیم مستقیماً از توزیع هدف \(f_X(x)\) نمونه‌گیری کنیم، از یک توزیع پیشنهادی ساده‌تر \(f_Y(x)\) کمک می‌گیریم.

قضیه ۴.۵.۴ (شرط پوشش)

باید یک مقدار ثابت \(c > 0\) وجود داشته باشد به طوری که برای تمام مقادیر \(x\):

\[f_X(x) \le c f_Y(x), \quad \forall x \in \mathbb{R}\]

این یعنی تابع چگالی پیشنهادی ضرب‌در \(c\) باید همواره بالاتر از تابع هدف باشد. بهترین \(c\) برابر است با ماکزیمم نسبت \(f_X(x)/f_Y(x)\).

الگوریتم شبیه‌سازی

  1. تولید کاندید: یک مقدار \(X^*\) از توزیع پیشنهادی \(Y\) تولید کنید.
  2. تولید یکنواخت: یک عدد \(U \sim U(0,1)\) تولید کنید.
  3. شرط رد/پذیرش: اگر شرط زیر برقرار بود، \(X^*\) را به عنوان نمونه نهایی بپذیرید:

\[U \le \frac{f_X(X^*)}{c f_Y(X^*)}\]

احتمال پذیرش برابر با \(\frac{1}{c}\) است. هرچه \(c\) به ۱ نزدیک‌تر باشد، الگوریتم بهینه‌تر است.

مثال: توزیع بتا (Beta(3,2))

در این مثال می‌خواهیم توزیع بتا را با استفاده از توزیع پیشنهادی یکنواخت شبیه‌سازی کنیم. تابع هدف: \(f_X(x) = 12x^2(1-x)\) برای \(0 \le x \le 1\). تابع پیشنهادی: \(f_Y(x) = 1\). ثابت بهینه: \(c = 16/9\).

کد زیر پیاده‌سازی دقیق این مثال است:

f_X <- function(x) {
  if (x >= 0 & x <= 1) {
    return(12 * (x^2) * (1 - x))
  } else {
    return(0)
  }
}

f_Y <- function(x) {
  if (x >= 0 & x <= 1) {
    return(1)
  } else {
    return(0)
  }
}

C <- 16 / 9

acbeta <- function(N) {
  samples <- numeric(N)  
  count <- 0             
  b <- 0 
  
  while (count < N) {
    b <- b + 1
    
    X_star <- runif(1, 0, 1)
    U <- runif(1, 0, 1)
    
    ratio <- f_X(X_star) / C
    
    if (U <= ratio) {
      count <- count + 1
      samples[count] <- X_star
    }
  }
  
 
  efficiency <- N / b

  
  N / b
  b
  round(efficiency, 4)
  round(1 / C, 4)
  
  return(samples)
}

acbeta(5)
## [1] 0.8060452 0.2664837 0.4223473 0.5749962 0.5866572
length(acbeta(5))
## [1] 5

پایان گزارش | روش‌های شبیه‌سازی آماری

```